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1 Report Summary 


1.1 Overview 

Reported here are our results of our numerical/theoretical investigation into 
the effects of thermal stress in nonisothermal gases under microgravity con- 
ditions. The first part of the report consists of a brief summary of the ac- 
complishments and conclusions of our work. The second part consists of two 
manuscripts, one being a paper presented at the 1998 MSAD Fluid Physics 
workshop, and the other to appear in Physics of Fluids. 

1.2 Project Objectives 

The objectives of our project were, first of all, to determine the accuracy of 
the Burnett constitutive equations for gases as applied to highly nonisother- 
mal, slow-moving gases by comparison ‘numerical experiments’ provided by 
direct simulation monte carlo (DSMC) methods. Secondly, we were to use 
these findings to assess the feasibility of using the microgravity environment 
to experimentally isolate and measure the flows resulting from thermal stress 
(as predicted by the Burnett equations) in highly nonisothermal, buovancy- 
free gases. 

2 Technical Summary of the Project 

2.1 Background and Motivation 

The continuum description of momentum and energy transport in gases, 
based upon Newton— Stokes— Fourier constitutive relations, can become in- 
accurate in rarefied or highly nonequilbrium regimes, i.e.. regimes in which 
the Knudsen number Kn (= X/L, where A is the gas mean free path and L 
is the characteristic system or gradient length) is no longer small. The Bur- 
nett equations, which represent the order— Kn 2 solution to the Boltzmann 
equation, ostensibly provide a means of extending continuum formulations 
into the transitional Knudsen regimes (~ Kn < 0.1). The accuracy and va- 
lidity of the Burnett equations, however, have not been firmly established. 
As has been noted by several authors, the asymptotic series expansion of 
the molecular distribution function — from which the Burnett equations are 
derived — has unknown convergence properties for finite Kn. The Burnett 
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equations can also lead to Second— law impossibilities, such as heat flux in 
an isothermal gas. Furthermore, the Burnett equations increase the order 
of the differential equations that govern momentum and heat transport in 
the gas. Additional boundary condition information is required to fully close 
the problem — yet such information is generally not available from physical 
principles alone. 

Because of these issues, it is generally held that the Burnett equations 
are valid only in regimes in which the Navier-Stokes-Fourier level of approxi- 
mation already provides an adequate description of transport, i.e. , regimes in 
which the Burnett contributions represent a small perturbation to heat and 
momentum transport. Such conditions can be representative of high— Mach 
number flows, for which application of the Burnett equations appears to 
have been the most successful. On the other hand, there is not a broad 
understanding of the accuracy of the Burnett equations when applied to 
slow— moving, nonisothermal flow conditions. In principle, ‘thermal stresses’ 
(fluid stresses resulting from temperature gradients — which are predicted 
by the Burnett equations) could become a significant convection mechanism 
in buoyancy— free, nonisothermal gases. Indeed, it has been recently sug- 
gested that thermal stress convection could affect the growth of crystals in 
microgravity physical vapor transport experiments. 

The project described here consisted of a theoretical and numerical ex- 
amination of thermally— induced stresses and flows in enclosed, highly non- 
isothermal gases under buoyancy— free conditions. A central objective has 
been to identify a strategy in which stress and/or convection effects, as pre- 
dicted by the Burnett equations, could be isolated and measured in micro- 
gravity— based experiments. Because of the questionable veracity of the Bur- 
nett equations, a second objective has been to test Burnett predictions of 
nonisothermal gas stress and convection with the exact description provided 
by the direct simulation Monte Carlo (DSMC) method. 

2.2 Thermal stress in two-dimensional gases 

The initial phase of the project was aimed at calculation, using continuum 
and DSMC methods, of gas convection in two dimensional nonuniformly 
heated rectangular enclosures. Typically, two adjacent surfaces of the en- 
closure were modeled as adiabatic, zero— stress surfaces (i.e., planes of sym- 
metry). and the other two adjacent surfaces were maintained at specified 
temperature distributions with one surface transferring a net amount of heat 
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to the gas, and the other transferring the heat from the gas. According 
to continuum theory, convection in the enclosure would result from thermal 
stress in the bulk gas and thermal creep along the side walls — the latter 
mechanism being the slip flow of gas over a solid surface that is driven by a 
gas temperature gradient tangential to the surface. 

Our continuum and DSMC calculations to date have not identified con- 
ditions in the enclosure that lead to measurable thermal stress flows that 
are comparable to or larger than thermal creep flows, and simultaneously 
maintain the Kn <0.1 regime required of the Burnett equations. With 
the exception of the pure continuum limit (Kn — > 0, under which thermal 
stress vanishes), elimination of thermal creep cannot be accomplished by 
maintaining the heated/cooled walls at uniform temperatures. Rather, the 
discontinuity (or jump) between the surface and adjacent gas temperatures — 
which will be proportional to Kn and the local normal temperature gradient 
— will lead to nonuniform gas temperatures along the nonuniformly heated 
surfaces. For all realistic values of Kn, thermal creep flows generated by the 
temperature jump effects were substantially larger than those resulting from 
thermal stress. 

To minimize the effects of creep, we performed additional simulations in 
which the temperature distributions along the heated/cooled surfaces were 
assigned to provide, for a given Kn, nearly uniform gas temperature adjacent 
to the surface. Surface temperature distributions were determined from so- 
lution of the gas conduction equation with uniform gas temperature bound- 
ary conditions along the adjacent heated/cooled surfaces, and subsequent 
application of the solution into the order— An temperature jump relations. 
This approach imposed a surface temperature variation along the heated wall 
which increased towards the junction with the cooled wall, with an opposite 
trend along the cooled wall. The effect of this strategy largely eliminated 
thermal creep from continuum— based models, and left a computed flow field 
which was driven almost entirely by thermal stress. However, this strategy 
would also lead to highly nonequilibrium conditions in the vicinity of the 
hot/cold surface junction, for which the Burnett equations are not expected 
to hold. DSMC calculations on the same system showed convective flows that 
were qualitatively similar yet substantially larger than those predicted by the 
continuum model — even though the DSMC and continuum— calculated tem- 
perature fields were in significant agreement. 
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2.3 Thermal stress in one-dimensional gases 

A more direct method of examining thermal stress effects would be to di- 
rectly measure the normal stress in a gas that is contained between two 
parallel surfaces at different yet uniform temperature — which will eliminate 
convection from thermal creep. We have investigated the effects of thermal 
stress, as predicted by the Burnett equation, on the pressure distributions 
and normal stress in a stationary, buoyancy— free, hard— sphere gas for the 
case of one— dimensional heat transfer. Using First— Law principles and the 
Burnett equation, it can be shown that thermal stress results in a reduction 
in normal stress in the nonisothermal gas relative to that in the equilibrium 
state. The normal stress, in turn, can be obtained as an eigenvalue to a 
second— order ordinary differential equation, representing the Burnett equa- 
tion, for the pressure distribution in the gas. Boundary conditions for the 
pressure at the heated and cooled surfaces can be obtained by an asymptotic 
expansion of the Burnett equations within the Knudsen rarefaction layers 
adjacent to the surfaces, which provides an order-Kh 2 correction the order- 
Kn pressure slip relations that are obtained from solution of the linearized 
Boltzmann equation at the gas/surface interface. 

The Burnett and DSMC predictions of pressure are in close agreement for 
effective Knudsen numbers (based on the temperature gradient in the gas) 
less than 0.1. In particular, the Burnett equations can accurately describe 
the shape of the Knudsen (or rarefaction) layers adjacent to the heated and 
cooled surfaces that bound the gas, and can also describe the variation in 
pressure in the ‘bulk’ gas (i.e., outside the Knudsen layers). In addition, 
theoretical predictions of the reduction in normal stress correspond well to 
DSMC— derived values. 

2.4 Conclusions 

2.4.1 Validity of the Burnett equations 

The original project was motivated by the concept of a slow-moving, non- 
isothermal gas flow driven entirely by thermal stress within the gas. Since 
buoyant forces would typically overwhelm the thermal stress forces in normal 
gravity conditions, we hypothesized, in the original proposal, that thermal 
stress flows could be realized under microgravity conditions. The central 
objective of the project was to conduct numerical simulations to determine 
1) the validity of the Burnett equations as applied to slow-moving thermal 
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stress flows, and 2) ascertain the feasibility of using the microgravitv envi- 
ronment to experimentally measure such flows. 

Regarding objective 1) above, we were successful in demonstrating that 
the Burnett equations are consistent with DSMC ‘reality’ for the simple case 
of a 1-D, stationary nonisothermal gas. On the other hand, we obtained 
decidedly ambiguous results on the validity of the Burnett equations applied 
to 2-D nonisothermal, slow-moving gases. For such situations, DSMC re- 
sults are only qualitatively similar to that predicted by the continuum-based 
Burnett model. 

A problem with the 2-D simulations - which we were not able to overcome 
- was that the creation of wall heating conditions necessary to drive thermal 
stress flows (as predicted by the Burnett equations) lead to local Knudsen 
conditions that were outside of the range of validity of the Burnett equations. 
In principle, this could be alleviated by reducing the overall Knudsen condi- 
tions of the simulation to the near-continuum level (i.e., Kn « 1) - yet it is 
currently not feasible to simulate such conditions with DSMC. Specifically, 
the near-continuum case would require a very large set of computational 
molecules (since the 2-D DSMC method scales as l/Kn ), and the very slow 
moving thermal stress flow in the near-continuum limit would be extremely 
difficult to resolve with DSMC statistics. 

2.4.2 Role of microgravity 

Microgravity conditions would not be needed to examine experimentally the 
effect of thermal stress on normal stress in 1-D, stationary nonisothermal 
gases. The hydrostatic pressure gradient would be insignificant under the 
Kn ~ 0.1 conditions necessary to resolve thermal stress effects. 

Because of the inability of current DSMC methods to model near-con- 
tinuum flows, we do continue to see a potential for using microgravity to 
examine the validity of the Burnett equations applied to 2-D nonisothermal 

gases. 


3 Project Accomplishments 

3.1 Graduate students 

Mr. Rong Wei (MS, 1999) performed DSMC and continuum calculations of 
hard-sphere gas heat transfer and fluid flow in two-dimensional, rectangular 
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enclosures. His thesis examined the extent to which the Navier-Stokes- 
Fourier formulation of transport, with slip and jump corrected boundary 
conditions, could model heat transfer in rarefied gases with two-dimensional 
temperature profiles. 

Mr. Joseph Ragan (MS) is currently working on the simulation of gas 
flow over heated surfaces that have directionally-dependent molecular ac- 
commodation properties. He was partially supported by the project for the 
first two quarters of 1999, and has recently received fellowship support from 
the National Defense Science and Engineering Grant program. 

3.2 Publications and Presentations 

1. “Drag and torque on iV-aggregated spheres in creeping flow,” A. V. 
Filippov, D. W. Mackowski, and D. E. Rosner, to be presented at the 
AIChE 1999 Annual Meeting. 

2. Mackowski, D. W., Papdopoulos, D. H., and Rosner, D. E., “Compari- 
son of Burnett and DSMC prediction of pressure distribution and nor- 
mal stress in one dimensional, strongly nonisothermal gases,” Physics 
of Fluids, to appear (1999). 

3. “Modelling of heat and mass transfer in rarefied, highly-nonequilibrium 
conditions,” Invited seminar, Mechanical Engineering Department, Uni- 
versity of Alabama at Birmingham. May 1999. 

4. “Investigation of thermal stress convection in nonisothermal gases un- 
der microgravity conditions,” presented at the 1998 MSAD Fluid Physics 
workshop, Cleveland, OH. 

5. “Modelling of heat and mass transfer in rarefied, highly-nonequilibrium 
conditions,” Invited seminar, Mechanical Engineering Department, Michi- 
gan State University, November 1997. 
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1. INTRODUCTION 

The continuum description of momentum and en- 
ergy transport in gases, based upon Newton -Stokes— 
Fourier constitutive relations, can become inaccurate in 
rarefied or highly nonequilbrium regimes, i.e., regimes 
in which the Knudsen number Kn (= A/L, where A is 
the gas mean free path and L is the characteristic system 
or gradient length) is no longer small. The Burnett 
equations, which represent the order— Kn 2 solution to 
the Boltzmann equation, ostensibly provide a means of 
extending continuum formulations into the transitional 
Knudsen regimes (~ Kn < 1 ). 

The accuracy and validity of the Burnett equations, 
however, have not been firmly established. As has 
been noted by several authors, the asymptotic series 
expansion of the molecular distribution function — from 
which the Burnett equations are derived — has unknown 
convergence properties for finite Kn , 1,2 The Burnett 
equations can also lead to Second — law impossibilities, 
such as heat flux in an isothermal gas . 3 Furthermore, the 
Burnett equations increase the order of the differential 
equations that govern momentum and heat transport in 
the gas. Additional boundary condition information is 
required to fully close the problem - yet such informa- 
tion is generally not available from physical principles 
alone. 

Because of these issues, it is generally held that 
the Burnett equations are valid only in regimes in which 
the Navier-Stokes-Fourier level of approximation al- 
ready provides an adequate description of transport, i.e., 
regimes in which the Burnett contributions represent 
a small perturbation to heat and momentum transport. 
Such conditions can be representative of high-Mach 
number flows, for which application of the Burnett 
equations appears to have been the most successful . 4-7 
On the other hand, there is not a broad understanding 
of the accuracy of the Burnett equations when applied 


cal and numerical examination of thermally-induced 
stresses and flows in enclosed, highly nonisothermai 
gases under buoyancy— free conditions. A central objec- 
tive has been to identify a strategy in which stress and/ or 
convection effects, as predicted by the Burnett equations, 
could be isolated and measured in microgravity— based 
experiments. Because of the questionable veracity of 
the Burnett equations, a second objective has been to 
test Burnett predictions of nonisothermai gas stress and 
convection with the exact description provided by the 
direct simulation Monte Carlo (DSMC) method. 

2. PREDICTION OF THERMAL STRESS 
CONVECTION 

The initial phase of the project was aimed at cal- 
culation, using continuum and DSMC methods, of gas 
convection in two dimensional nonuniformly heated 
rectangular enclosures. Typically, two adjacent surfaces 
of the enclosure were modeled as adiabatic, zero— stress 
surfaces (i.e., planes of symmetry), and the other two 
adjacent surfaces were maintained at specified temper- 
ature distributions with one surface transferring a net 
amount of heat to the gas, and the other transferring the 
heat from the gas. 

The continuum formulations of momentum and en- 
ergy transport are identical to Navier— Stokes— Fourier 
models, with the exception of the Burnett stress tensor 
in the momentum equations and the creep and jump 
boundary conditions. For the conditions examined here 
(i.e., slow— moving flow, with Rei <§C 1 ) the only sig- 
nificant terms in the Burnett stress tensor relations will 
be those involving temperature gradients. This thermal 
stress component appears as 8,11 

t t = -^ u , 3 (wr-i(V 2 T)f) 

+ |( ( VT ) (VT)-i(VT. v T) 1 ) (1) 


to slow— moving, nonisothermai flow (SNIF) condi- 
tions. As noted by Kogan, ‘thermal stresses’ (fluid 
stresses resulting from temperature gradients — which 
are predicted by the Burnett equations) could become 
a significant convection mechanism in buoyancy— free, 
nonisothermai gases . 8,9 Indeed, it has been recently 


in which /i is the dynamic viscosity, R is the gas 
constant, and CJ 3 and CU 5 are dimensionless, order— unity 
coefficients which depend on the interaction potential of 
the molecules. The creep and jump boundary conditions 
appear 


suggested that thermal stress convection could affect 
the growth of crystals in microgravity physical vapor 
transport experiments . 10,11 

The work presented here consists of a theoreti- T — T w -\ — - 7 ^ 


. dT\ 
U On) 


(3) 
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where n is the outward normal and thedimensionless co- 
efficients c s and ct depend on the thermal and momen- 
tum accommodation properties of the surface. 12 Nu- 
merical solution of the governing equations was accom- 
plished using the SIMPLER algorithm of Patankar. 13 
Coefficients corresponding to hard-sphere molecules, 
which gives a temperature-dependent viscosity of fi ~ 
T l/2 , were used in the computations. 

Direct simulation Monte Carlo calculations of 
hard-sphere gas convection and heat transfer were 
accomplished using the standard procedure developed 
by Bird. 14,15 The cell size was nominally set to 0. 1 A 0 , 
where A 0 represents the mean-ffee-path at the equi- 
librium state of the system, and 10—20 molecules were 
assigned per cell. Simulations were conducted for a 
Knudsen range of Kn = 0.01 - 0.2. Because thermal 
creep and stress flows will be on the order of Kn times 
the mean molecular velocity, resolution of the flows 
using DSMC required simulation times on the order of 
10 6 — 10 7 time steps. 

Our continuum and DSMC calculations to date in- 
dicate that it would be very difficult to create conditions 
in the enclosure that result in measurable thermal stress 
flows that are comparable to or larger than thermal 
creep flows, and simultaneously maintain the Kn < 0.1 
regime required of the Burnett equations. With the ex- 
ception of the pure continuum limit (Kn -> 0, under 
which thermal stress vanishes), elimination of ther- 
mal creep cannot be accomplished by maintaining the 
heated/ cooled walls at uniform temperatures. Rather, 
the discontinuity (or jump) between the surface and ad- 



jacent gas temperatures - which will be proportional to 
Kn and the local normal temperature gradient - will 
lead to nonuniform gas temperatures along the nonuni- 
formly heated surfaces. For all realistic values of Kn, 
thermal creep flows generated by the temperature jump 
effects were substantially larger than those resulting 
from thermal stress. 

To minimize the effects of creep, we performed 
additional simulations in which the temperature dis- 
tributions along the heated/ cooled surfaces were as- 
signed to provide, for a given Kn y nearly uniform gas 
temperature adjacent to the surface. Surface tempera- 
ture distributions were determined from solution of the 
gas conduction equation with uniform gas temperature 
boundary conditions along the heated/ cooled surfaces, 
and subsequent application of the solution into Eq. (3) 
to predict T w (x), This approach imposed a surface 
temperature on the heated surface which increased to- 
wards the junction with the cooled surface, with an 
opposite trend along the cooled wall. The effect of 
this strategy resulted in thermal creep flows that were 
confined about the hot/cold junction, and left a bulk, 


thermal-stress driven convection pattern in the bulk 
of the enclosure. However, this strategy would also 
lead to highly nonequilibrium conditions in the vicinity 
of the hot/ cold surface junction, for which the Burnett 
equations are not expected to hold. 

The comparison between continuum/ Burnett and 
DSMC predictions of convective flows in the enclo- 
sure has been inconclusive. DSMC calculations have 
shown convective flows that are qualitatively similar 
than those predicted by the continuum model. As an 
example, we show in Fig. 1 plots of velocity vec- 
tors and isotherms, calculated using the DSMC (top) 
and continuum (bottom) models, for a Kn = 0.02 
hard-sphere gas contained in an enclosure with a rela- 
tive temperature difference on the hot and cold walls of 
2 (Th — Tq)/(Th -f Tc) = 1. The heated wall is on the 
top, and the left and bottom walls are symmetry surfaces. 
A strip of the walls adjacent to the hot/ cold junction, 
equal to 0. 1 of the wall length, is held adiabatic, and the 
remainder of the walls have temperature distributions 
set to give isothermal gas conditions per the procedure 
discussed above. Velocity vectors corresponding to ther- 
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mai creep along the adiabatic surfaces (seen in the tight 
counterclockwise rotation in the upper right corner) have 
been removed from the plot to allow resolution of the 
thermal stress flows. The stress flow - as predicted by 
the continuum model - results in a clockwise rotation 
in the main body of the gas for the given conditions. 
A similar pattern is seen in the DSMC results. We 
cannot establish, however, whether the observed DSMC 
flows result from thermal stress, or are due to slip ef- 
fects at the walls. As mentioned above, the veracity of 
the Burnett relations is also questionable for the highly 
nonequilibrium conditions of the simulations. On the 
other hand, temperature profiles calculated via DSMC 
and continuum models show significant agreement. 

3. THERMAL STRESS IN 1-D HEAT 
TRANSFER 

A simpler situation in which to compare contin- 
uum/Burnett and DSMC predictions of thermal stress 
effects is offered by 1-D heat transfer in a stationary 
gas. In this situation, the effects of thermal stress would 
be seen in the pressure distribution and normal stress in 
the gas. 16 

The computational domain was now taken to be 
a slab of gas contained between two parallel surfaces, 
separated by a distance L, with the surfaces at x = 0 


temperature will be described by 

q* = constant = 8 1 ^ 2 0 t (5) 

where q * is the dimensionless heat flux (= qL/kT m ). 
Equation (5) can be used to combine the last two terms 
in Eq. (4), which results in 



where c<i — — 2cj 5 )/2 = 0.9900 for hard-sphere 

molecules. 

Two separate effects — or regimes - on 0 can be 
anticipated from inspection ofEq. (6). One effect, which 
is discussed by Kogan 18 and Makashev 19 , derives from 
the fact that the derivatives of 0 will be multiplied by 
the small parameter (for near-continuum conditions) of 
Kn 2 . The solution to Eq. (6) could therefore exhibit 
‘boundary layers’ of width ~ Kn. It is shown below 
that this property, combined with appropriate boundary 
conditions, will allow for a limited description of the 
Knudsen layers adjacent to the surfaces. 


and L maintained at uniform temperatures of Tcs and 
Ths (with Ths > Tcs)> respectively. In nondimen- 
sional form (with pressure and stress normalized with 
the equilibrium pressure P m and temperature by the 
equilibrium temperature T m ), the Burnett equation for 
the x-directed, normal component of the stress tensor 


* , C\Kn 2 6 

T = 0 + — 

0 


9*i >' 

0J4 — CJ2 

0 




A second characteristic regime, as indicated by 
Eq. (6), would occur outside the Knudsen layers. By 
expanding 0 in a power series of Kn 2 , and neglecting 
all terms higher than Kn 2 , the pressure distribution in 
the bulk gas would be given approximately by 


0 % r* 


C\C 2 (Kn q*) 2 
~8 


(7) 


As is evident from inspection of Eq. (7), thermal stress 


0/2 ' 

+u)$Q” + u)$— ( 4 ) 

In the above, 0 = P/P m , r* = r/P m , 9 = T/T m , the 
prime denotes differentiation with respect to f = x/L, 
C\ = (47t/ 3)(5/16) 2 = 0.4091, and the dimensionless 
cj coefficients depend on the molecular interaction po- 
tential. Since the gas is stationary and buoyancy— free, 
the stress r will be a constant. In the limit of Kn 0, 
this gives the Navier— Stokes result of P — r = con- 
stant. For finite Kn , however, the additional source 
of thermal stress can act within the nonisothermal gas. 
The magnitude of the thermal stress will vary with po- 
sition — by virtue of the dependence of temperature and 
temperature gradient on position - and consequently 
pressure will vary to maintain a constant normal stress. 

The Burnett equations make no contribution to the 
heat flux for a stationary gas. Consequently, the gas 


would create a pressure gradient in the gas, with pres- 
sure increasing towards the cooler regions in the gas. 
The gradient would be proportional to the square of 
Knq* ~ \/8d6/d(xf A) — which can be interpreted as 
a Knudsen number based on the characteristic length of 
the temperature gradient (note that this quantity is inde- 
pendent of L). It should be emphasized that the effect 
predicted from Eq. (7) is fundamentally different than 
the pressure gradient created by ‘thermal transpiration’ 
of a gas in a tube with an imposed axial temperature 
gradient. 20,21 The latter is a result of thermal slip at the 
walls ofthe tube, and leads to a pressure that increases in 
the direction of increasing temperature. Thermal stress, 
on the other hand, results from the effect of temperature 
gradients on the molecular velocity distribution function 
within the gas. 

Although the thermal stress pressure gradient can 
be labeled ‘hydrostatic’ - since the gas is at rest - it 
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is distinctly different than that resulting from a gravita- 
tional acceleration in the x— direction. Unlike the latter, 
thermal stress would not result in a difference between 
the normal forces acting on the hot and cold surfaces. In 
other words, the ‘pressure' measured at the surfaces — 
which would physically represent the normal stress r — 
would be identical for both surfaces. 

Thermal stress, however, will result in a different 
value of the normal stress than that predicted from the 
Navier-Stokes level of approximation. This follows 
from conservation ofenergy requirements. In particular, 
the average pressure in the gas represents the equilib- 
rium pressure that would be attained if the walls were 
instantaneously made adiabatic. Since the equilibrium 
pressure is used to normalize the dimensional pressure 
P, this statement is equivalent to 



Regardless of the values of q* and Kn, the pressure 
distribution in the gas must satisfy the energy conser- 
vation constraint implied by Eq. (8). Consequently, the 
normal stress r* would be obtained as the eigenvalue 
to Eq. (6) such that the solution (for specified boundary 
conditions) satisfies Eq. (8). In general, this value will 
be different than the Navier-Stokes result of r* = 1. 

An approximate value for r* can be obtained by 
neglecting the effects of the Knudsen layers at the 
surfaces, for which the pressure distribution would be 
given by Eq. (7). To order Kn 2 , this gives 16 

r* se 1 - cic 2 {Knq') 2 (9) 

This relatively— simple approximation indicates that 
thermal stress will lower the normal stress in a closed 
system relative to that predicted from the Navier-Stokes 
level - although we note again that the effects of Knud- 
sen layers have been neglected in the analysis. 

The final elements required to close the prob- 
lem are the boundary conditions for pressure. As is 
the case with the Navier-Stokes approximation, the 
boundary conditions for the Burnett equations should 
represent an extrapolation of the solution across the 
region, adjacent to the wall, where the solution is no 
longer valid. Makashev 19 and Schamberg 22 have pro- 
posed boundary conditions that are consistent with the 
order— /f n 2 accuracy of the Burnett equations. The 
accuracy of these approaches, however, has not been 
well established. 23 Alternatively, order -Kn relations 
can be derived for the pressure ‘slip’ adjacent to a heated 
or cooled surface. 12,21 However, our work at this stage 
is primarily concerned with determining whether there 
are boundary conditions which, when coupled to the 



Figs. 2-4: DSMC and continuum pressure distributions 

Burnett equations, can reproduce DSMC predictions of 
pressuredistributions in a nonisothermal gas. Therefore, 
the pressures at the hot and cold surfaces were taken 
to be parameters, and were chosen to provide the best 
agreement between theory and DSMC results. The ob- 
vious choice for the pressure at the surfaces will be the 
values determined from DSMC predictions. 

Comparisons of Burnett (via numerical solution of 
Eq. (6)) and DSMC predictions of pressure distribution 
appear in Figs. 2-4. Each plot shows P/r = <j)/r* vs. 
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xi — xjL for a fixed value of q* , with Kn a parame- 
ter. Two theoretical predictions of P/r are shown for 
each set of DSMC results. The first corresponds to the 
numerical solution of Eq. (6), with boundary conditions 
obtained from exi rapolat ion of t he DSMC —derived pres- 
sures to the surfaces and temperatures predicted from 
Eq. (5). The second represents the bulk gas thermal stress 
pressure distribution predicted from Eq. (7), which does 
not account for boundary effects. This latter predic- 
tion has been shifted by a constant to match with the 
full-Burnett solution at f = 0.5. 

As is evident from the results, the pressure profiles 
show distinct Knudsen layers at both the cooled and 
heated surfaces. The drop in pressure at the cooled 
surface, and the increase in pressure at the hot surface, 
are both consistent with the predictions of pressure slip 
relations. 21 The pressure drop at the cold surface can 
be considerable for the conditions examined here — 
amounting to around an 8% decrease for Kn = 0.2 and 
q* - 1.5. 

The solutions of the Burnett equation, with DSMC 
derived boundary conditions, are seen to capture the 
essential features of the DSMC pressure distribution. 
In particular, the solutions provide a good description 
of the width and form of the Knudsen layers and the 
pressure distribution outside the layers. The difference 
between the theoretical and DSMC results is greatest at 
the edge of the cold-surface Knudsen layers, for which 
the theoretical model tends to overpredict the pressure. 
This is most evident for the results corresponding to 
q* = 2.0 in Fig. 4. Nevertheless, the fact that the 
Burnett equations can resolve, to a reasonable accuracy, 
the Knudsen layers at the surface is somewhat surprising 
- especially when considering that the theory is based on 
the order — Kn continuum temperature profile. We also 
examined solut ions to Eq. (6) using boundary values of <p 
that were different than the DSMC results, and found that 
the exact, DSMC— derived boundary conditions provide 
the best overall agreement between Burnett equation 
predictions and DSMC results. 

The DSMC results for q* = 2.0 appear to show 
a pressure distribution in the bulk gas that is described 
by Eq. (7). On the other hand, the pressure distribution 
for q * = 1.5 and Kn = 0.2 (Fig. 3) — for which 
Eq. (7) predicts a greater effect - is dominated by 
the Knudsen layers extending from the surfaces. To 
eliminate the effects of the Knudsen layers at the hot 
wall, we performed additional DSMC calculations in 
which the velocities of the incoming molecules at the 
hot boundary were sampled from the Chapman-Enskog 
distribution function for the fixed values of q* and Kn. 

By doing so, the hot surface now approximated an open 
boundary. The corresponding DSMC results showed a 



pressure distribution in the bulk gas that was accurately 
represented by Eq. (7). 

A final comparison of theory and experiment can 
be obtained from the dimensionless normal stress. The 
simplified model of Eq. (9) indicates that t* should be 
a function primarily of Knq* . Accordingly, we plot in 
Fig. 5 the DSMC values of r* vs. Knq* for the seven 
different combinations of Kn and q* that were used 
in the closed system calculations (results of Figs. 2-4). 
Theoretical results correspond to the derived eigenvalues 
of Eq. (6) for the DSMC-derived boundary conditions, 
and to the approximation given by Eq. (9). 

The first point to make is that the predictions of r* 
from full solution ofthe Burnett differential equation are 
nearly equivalent to those obtained from the bulk-gas 
approximation of Eq. (9). Evidently, the decrease in 
pressure at the cold surface is compensated by the 
increase at the hot, so that the Knudsen layers have 
a small effect on the averaged pressure in the gas. 
Secondly, the primary dependence of r* on Knq* is 
supported in the DSMC results at Knq* = 0.1 and 0.2 
- which each correspond to two combinations of Kn 
and q*. As observed, the results are nearly identical 
at these points. Finally, the theoretical predictions 
are in excellent agreement with the DSMC results for 
Knq* < 0.15, beyond which the theory overpredicts 
the decrease in r\ As can be seen from the results, the 
relative decrease in normal stress on the surfaces is quite 
small, i.e., r* = 0.975 for q* = 2.0 and Kn = 0.2, or 
a 2.5% decrease in ‘measured’ pressure at the surface. 
We should emphasize, however, that this decrease is 
still significantly larger than the numerical precision of 
the DSMC simulations. 

SUMMARY 

The project has sought to ascertain the veracity of 
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the Burnett relations, as applied to slow moving, highly 9 
nonisothermal gases, by comparison of convection and 
stress predictions with those generated by the DSMCio 
method. The Burnett equations were found to provide 
reasonable descriptions of the pressure distribution and 
normal stress in stationary gases with a I -D temperature 
gradient. Continuum/ Burnett predictions of thermal 
stress convection in 2— D heated enclosures, however, 
are not quantitatively supported by DSMC results. For 
such situations, it appears that thermal creep flows, 
generated at the boundaries of the enclosure, will be 12 
significantly larger than the flows resulting from thermal 
stress in the gas. 13 
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The Burnett equations have been shown to provide improved descriptions, relative 
to the Navier-Stokes equations, of flow structure in high-velocity (i.e., hypersonic) gases. 
We examine here the accuracy of the Burnett constitutive equation for fluid stress as ap- 
plied to stationary gases. Specifically, we investigate the effects of ‘thermal stress’ (fluid 
stress induced by a temperature gradient), as predicted by the Burnett equation, on the 
pressure distributions and normal stress in a stationary, buoyancv-free, hard-sphere gas 
for the case of one-dimensional heat transfer. We show, using First-Law principles and 
the Burnett equation, that thermal stress results in a reduction in normal stress in the 
nonisothermal gas relative to that in the equilibrium state. The normal stress, in turn, 
can be obtained as an eigenvalue to a second-order ordinary differential equation, repre- 
senting the Burnett equation, for the pressure distribution in the gas. Simple asymptotic 
solutions to the Burnett equation are developed, and are used in combination with order- 
Kn pressure slip relations to formulate pressure boundary conditions at the heated and 
cooled surfaces. The approximate solutions, as well as exact numerical calculations, are 
compared with pressure distributions generated from the direct-simulation-Monte-Carlo 
(DSMC) method. The Burnett and DSMC predictions of pressure are in good agreement 
for effective Knudsen numbers (based on the temperature gradient in the gas) less than 
0.1. In particular, the Burnett equations can provide a reasonable description of the Knud- 
sen (or rarefaction) layers adjacent to the heated and cooled surfaces that bound the gas, 
and can also describe the variation in pressure in the bulk gas. In addition, theoretical 
predictions of the reduction in normal stress correspond well to DSMC-derived values. 
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I. BACKGROUND AND MOTIVATION 

In 1938. Burnett developed an approximate solution to the Boltzmann equation, based 
upon the asymptotic expansion procedure of Chapman and Enskog, that provided consti- 
tutive relations accurate to order Kn 2 (where Kn = X/L is the Knudsen number, with 
A and L the mean free path and characteristic system length). 1 ’ 2 Ostensibly, these re- 
lations, commonly referred to as the Burnett equations, offered a means of extending 
continuum-based descriptions of momentum and heat transport in a gas to transitional 
Knudsen regimes. However, the accuracy and validity of the Burnett equations have not 
been firmly established. As has been noted by several authors, the asymptotic series ex- 
pansion of the molecular distribution function - from which the Burnett equations are 
derived - has unknown convergence properties for finite Kn. 2 ’ 3 Furthermore, the Burnett 
equations increase the order of the differential equations that govern momentum and heat 
transport in the gas. Additional boundary condition information is required to fully close 
the problem - yet such information is generally not available from physical principles alone. 
Lee has recently shown that, as expected, the lack of higher— order boundary conditions, 
coupled with the Burnett equations, can result in non-unique solutions of velocity, tem- 
perature, and pressure in steady Couette flow. The Burnett equations can also lead to 
Second-Law impossibilities in certain situations, such as a negative dissipation function or 
a heat flux in an isothermal gas. 5,6 

Because of these issues, it is generally held that the Burnett equations are valid only 
in regimes in which the Navier-Stokes-Fourier level of approximation already provides an 
adequate description of transport, i.e., regimes in which the Burnett contributions rep- 
resent a small perturbation to heat and momentum transport. Such conditions can be 
representative of high-Mach number flows, for which application of the Burnett equations 
appears to have been the most successful. 7-10 On the other hand, there is not a broad 
understanding of the accuracy of the Burnett equations when applied to slow-moving, non- 
isothermal flow conditions. As noted by Kogan, ‘thermal stresses’ (fluid stresses resulting 
from temperature gradients - which are predicted by the Burnett equations) could become 
a significant convection mechanism in buoyancy-free, nonisothermal gases. 11,12 Indeed, 
it has been recently recognized that thermal stress convection could affect the growth of 
crystals in highly-nonisothermal yet buoyancy-free microgravity physical vapor transport 
experiments. 13,14 

An estimation of the characteristic thermal stress velocity in slow-moving, nonisother- 
mal flow conditions can be obtained from an order-of-magnitude analysis of the Burnett 
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equations, which yields 14 


AT 

~T~ 
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1 / 

Vts ~ ^ 


where v is the gas kinematic viscosity, and T and AT are the characteristic system temper- 
ature and temperature difference, respectively. Assuming (A T)/T ~ 1, the thermal stress 
velocity will be on the order of the momentum diffusion velocity - which corresponds to 
order-unity Reynolds flows. Such ‘creeping’ flows would, in principle, represent a small 
perturbation from the zero velocity and constant pressure conditions that would be triv- 
ially predicted by the Navier-Stokes equations (with no-slip boundary conditions). Yet 
by this very reasoning, the presence of flows and pressure gradients in a nonisothermal, 
buoyancy-free gas would be an unambiguous indication of the effects of thermal stress. 

The objective of this paper is to, first of all, apply the Burnett equations to the 
situation of one-dimensional heat transfer in a quiescent gas, and theoretically predict 
the pressure distribution that is induced in the gas from thermal stress. Secondly, we will 
employ the ‘numerical experiment’ provided by the direct simulation Monte Carlo (DSMC) 
method, applied to the same system, to ascertain the accuracy of the Burnett equation 
predictions. We focus on pressure effects in a 1-D, stationary, nonisothermal gas in the 
absence of gravity, as opposed to 2-D gas motion induced by thermal stress, primarily 
because the weak, order-An 2 convection predicted by the Burnett equations would be 
difficult to resolve with DSMC statistics. In addition, a 2-D nonisothermal configuration 
would result in gas convection from thermal slip at the boundaries, which would interfere 
with our intention of isolating the effects of thermal stress. 15,16 

Investigations along the same lines as those conducted here have been recently per- 
formed by Wadsworth 17 and Zhong and Koura 18 , who used DSMC to examine the accuracy 
of continuum formulations (Navier-Stokes with or without the Burnett terms) for 1-D heat 
transfer. In all previous examinations of which we are aware, the ability of the Burnett 
equations to predict pressure distributions in a quiescent, non-isothermal gas has not been 
addressed. One reason why this effect has escaped scrutiny, as will be shown below, is 
that it is essentially negligible under moderately nonisothermal and near-continuum (e.g., 
An < 0.1) conditions. A second reason has been the lack of well-established boundary 
conditions for pressure at a heated/cooled surface (i.e., pressure slip) which are consistent 
with the order- Ah 2 accuracy of the Burnett equations. Our approach to this latter point 
will be to develop, via a simple asymptotic analysis of the Burnett stress equation in the 
Knudsen layer, a thermal stress correction to order-An pressure slip at the surfaces. 
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II. THEORETICAL FORMULATION 


A. System Configuration and Problem Definition 

The system, which is schematically illustrated in Fig. 1, consists of a stationary, hard 
sphere monatomic gas that is contained between two parallel, infinite-area boundaries that 
are separated by a distance L. A heat flux of q is flowing in the gas in the negative x 
direction. The system is in steady-state, and gravity is absent. 

To identify the relevant system parameters which affect the pressure distribution, we 
present our analysis in a nondimensional form. The nondimensional temperature, pressure, 
and position are defined 



0 = 


P 

K 7 f ' 



(i) 


where T re f and P re f are an appropriate reference temperature and pressure, which are 
defined below. We use the viscosity-based mean free path to define the Knudsen number; 

1/2 
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where R is the specific gas constant and n is the dynamic viscosity. The thermal con- 
ductivity k is given by k — 15f?/i/4, and both k and /i are proportional to T 1 / 2 for the 
hard-sphere gas. Finally, the nondimensional heat flux and normal fluid stress are defined 
as 


Q_L 

Kef T re f 



(3) 


B. Continuum Formulation 


In nondimensional form, the Burnett equation for the a: -directed, normal component 
of the stress tensor is 4,10 


= </> + 
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where the prime denotes differentiation with respect to £, c\ = (47t/ 3)(5/16) 2 = 0.4091, 
and the dimensionless ui coefficients depend on the molecular interaction potential. For 
hard-sphere molecules, the values are 


u 2 = 2.028, w 3 = 2.418, u 4 = 0.681. u; 5 = 0.219 (5) 

We include the hydrostatic pressure P into our definition of the normal stress solely for 
convenience. Since the gas is stationary and buoyancy-free, the stress r will be a constant. 
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In the limit of Kn — > 0, this gives the Navier-Stokes result of P — r — constant. For finite 
Kn, however, the additional source of thermal stress can act within the nonisothermal gas. 
The magnitude of the thermal stress will vary with position by virtue of the dependence 
of temperature and temperature gradient on position - and consequently pressure will vary 
to maintain a constant normal stress. 

The Burnett equations make no contribution to the heat flux for a stationary gas. 
Consequently, the gas temperature will be described by 

q* = constant = 9 l ^ 2 9' (6) 


This is integrated from £ = 0 (the cold surface) to yield 


9 = 




( 7 ) 


As is well recognized, the dimensionless temperature 9c should not be interpreted as the ac- 
tual gas temperature directly at the surface. Rather, it represents the extrapolation to the 
surface of the temperature distribution outside the Knudsen (or rarefaction) layer. 17,19 20 
This ‘apparent’ gas surface temperature is related to the solid surface temperature, 9cs, 
by a slip condition, which appears as 


9c — ®cs (l + c t Kn9 £ -+o) 


( 8 ) 


in which c t is a constant. We use the value c t = 1.691, which was derived by Loyalka by 
solution of the linearized Boltzmann equation for 1— D heat transfer in a hard-sphere gas 
and for perfectly accommodating surfaces. 19 By setting £ = 1 in Eq. (7), the heat flux is 
related to the hot-side boundary temperature by 

= | ($ 2 - 2 ) 

The corresponding slip condition at the hot surface is 

9 h = 9hs (l ~ c t Kn9 l £_ yl ) (10) 

For reasons which will become apparent in a following section, we will now define the 
reference temperature T re / so that 
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#c= 24 


By combining Eqs. (9) and (11), #c and @H can be obtained as functions of q , i.e., 

24 + q* 2 - \/3 q* (48 - q * 2 ) ^ ( 12 ) 

Note that the maximum value of q* , which corresponds to 9c -> 0, is 2\/3 = 3.464. 

Equation (6) can be used to combine the last two terms in Eq. (4), which results in 

ciKn 2 \ al uy \ _ w | ( 13 ) 


r* = 0 + 


0 


'( W (W\'\ <*** 


where c 2 = (w 3 - 2w 5 )/2 = 0.9900. Equation (13) provides a second-order, nonlinear 
ordinary differential equation for the pressure distribution in the gas. Complete solution 
requires specification of the nondimensional stress t* and a pair of boundary conditions for 
0. Alternatively, the equation could be differentiated, which would eliminate r* yet add 
an additional boundary condition for 0 - since the equation would become third order. 
We will work with the equation as it appears in Eq. (13), and employ thermodynamic 
principles to infer the value of t* . Our approach is discussed in the next section. 


C. Characteristics of Thermal Stress and Boundary Conditions 


Equation (13) displays the characteristic stiff form of the Burnett equations, in which 
the derivatives of the dependent variables (in this case 0) are multiplied by the small 
parameter An 2 . Accordingly, the solution to Eq. (13) would be expected to display an 
inner (or boundary layer) regime, in which the derivatives of 0 contribute to the solution, 
and an outer regime in which the derivatives have little effect. Kogan 21 and Makashev 22 
have noted that this property of the Burnett equations, when combined with order -An 2 
boundary conditions, would allow for a limited description of the Knudsen layers that 
occur at the surface/gas interfaces. 

We will ultimately use a numerical solution to Eq. (13) to make quantitative com- 
parisons between the Burnett and DSMC predictions of pressure distribution. However, 
insight into the effects of thermal stress on the pressure and normal stress in the gas - as 
well as a check of the consistency of the solution with regard to An constraints - can be 
obtained from relatively simple, order-An 2 asymptotic approximations to the inner and 


outer solution regimes to Eq. (13). 

The outer regime is examined first. By expanding 0 in a power series of An 2 , and 
neglecting all terms higher than An 2 , the pressure distribution in the bulk gas would be 
given approximately by 

CiC 2 {Knq *) 2 


0 ~ T* ”b 
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(14) 
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As is evident from inspection of Eq. (14), thermal stress would create a pressure gradient 
in the gas, with pressure increasing towards the cooler regions in the gas. The gradient 
would be proportional to the square of Knq* ~ \/dd9/d(x/X) - which can be interpreted 
as a Knudsen number based on the characteristic length of the temperature gradient (note 
that this quantity is independent of L). Accordingly, the outer solution would be valid 
only for Knq* C 1. It should be emphasized that the effect predicted from Eq. (14) is 
fundamentally different than the pressure gradient created by ‘thermal transpiration’ of a 
gas in a tube with an imposed axial temperature gradient. 17 23 The latter is a result of 
thermal slip at the walls of the tube, and leads to a pressure that increases in the direction 
of increasing temperature. Thermal stress, on the other hand, results from the effect of 
temperature gradients on the molecular velocity distribution function within the gas. 

Even though the thermal stress pressure gradient can be labeled ‘hydrostatic’ - since 
the gas is at rest - it is distinctly different than that resulting from a gravitational acceler- 
ation in the rc -direction. Unlike the latter, thermal stress would not result in a difference 
between the normal forces acting on the hot and cold surfaces. In other words, the ‘pres- 
sure’ measured at the surfaces - which would physically represent the normal stress r 
would be identical for both surfaces. 

Thermal stress, however, will result in a different value of the normal stress than that 
predicted from the Navier -Stokes level of approximation. This follows from conservation 
of energy requirements. In particular, the average pressure in the gas, given by 



represents the equilibrium pressure that would be attained if the walls were instantaneously 
made adiabatic. If the reference pressure is defined as P re f = P m , the above equation is 
equivalent to 

l 

4>d(, = 1 (15) 

Regardless of the values of q* and Kn , the pressure distribution in the gas must satisfy 
the energy conservation constraint implied by Eq. (15). Consequently, the normal stress 
t* would be obtained as the eigenvalue to Eq. (13) such that the solution (for specified 
boundary conditions) satisfies Eq. (15). In general, this value will be different than the 
Navier-Stokes result of r* = 1. 

An approximate value for r* can be obtained by neglecting the effects of the Knudsen 
layers at the surfaces, for which the pressure distribution would be given by the outer 
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solution of Eq. (14). Integration of this equation over £ gives 
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where we have used the definition of the reference temperature per Eq. (11). 
Eq. (16) for t* yields 
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1 

2 
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Solving 


(17) 


which, to order Kn 2 , is 


r* ss 1 — cic 2 (Knq*) 2 


(18) 


This relatively-simple approximation indicates that thermal stress will lower the normal 
stress in a closed system relative to that predicted from the Navier-Stokes level - although 
we note again that the effects of Knudsen layers have been neglected in the analysis. 

We turn now to an examination of the Knudsen layer region adjacent to the cold 
surface (an analysis at the hot surface would follow an analogous approach). Our objective 
here is to obtain a simple description of the boundary layer behavior predicted by Eq. (13). 
To do this, we assume that the Knudsen layer introduces an order Kn perturbation in the 
pressure field - which is based the fact that pressure slip at a heated/cooled surface occurs 
at the order-Kh level . 17,19 Accordingly, the pressure distribution in the Knudsen layer is 
estimated by superimposing the average bulk effect, given by Eq. (14), with an order - Kn 
undetermined function / lt via 


(pinner = t* + Kn fi + Kn 2 b (19) 

in which b ~ c.^q* 2 / (t* 9c ), with 9c representing the average temperature in the Knudsen 
layer, i.e., the average from £ = 0 to 9cKn. We now replace this form into Eq. (13), 
substitute £ with the stretched variable rj = £/Kn, expand the equation in powers of Kn , 
and retain terms to order-Kh. This gives the simple differential equation 

d 2 f i 

= ° (20) 

where a = t * 2 / (ciu^c) is an order-unity constant. The surface boundary condition for 
pressure is assumed to be in a similar form as Eq. (19); 

0(0) = (pc — + Kn fic + Kn 2 b (21) 
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and the outer boundary condition has /j — » 0 for r/ — > oc. Upon solution of Eqs. (20) and 
(21), the inner pressure distribution approximation is 


dinner = T* + Kn 2 b + (0c ~ T* - Kl % 2 ft) exp [— y/aTj] 


( 22 ) 


A composite solution, which approximates the pressure distribution in both the inner and 
outer regimes, is obtained by combining Eqs. (22) and (14), and adjusting the result by a 
constant to maintain the 0(0) = 0c condition. This results in 


<P 


r* + An 2 b 




+ (0c - r * 



(23) 


Since 1 — 9c /@c ~ O(Kn), the outer part of this approximation is consistent, to order 
Kn 2 , to that in Eq. (14). 

The boundary layer effect is observed, in this simplified analysis, in the exponential 
decay of pressure from the boundary to the freestream value. The characteristic width of 
the layer will be A£ « Kn/yfa ~ KnQc, he., proportional to the local mean free path, as 
expected. In addition, the pressure gradient at the surface, which is approximated by 


d<P 

dt 




Wc-t*) ^ +0(Kn 2 ) 


is order unity since 0c — r* is order-An. Consequently, the solution maintains the con- 
straint that Kn 0' ^ 1 within the Knudsen layer, providing that Kn <C 1 . 

The final elements required to close the problem - and allow for a quantitative compar- 
ison between the Burnett equation and DSMC predictions - are the boundary conditions 
for pressure. In formulating the problem for the inner solution, we posed the bound- 
ary conditions simply as a specified pressure directly at the surface. More precisely, the 
boundary conditions should represent an extrapolation of the continuum solution across 
the A£ ~ Kn 2 region adjacent to the wall (the ‘Burnett sublayer) within which the Burnett 
equation is no longer valid. Makashev 22 and Schamberg 24 have derived boundary condi- 
tions that are consistent with the order-An 2 level of the Burnett equations. However, the 
accuracy of these approaches has not been rigorously established 18 , nor are we aware of 
order-An 2 relations specifically for the pressure ‘slip’ at a heated or cooled surface. It 
should be noted that in application of the Burnett equations to 1-D hypersonic flow - 
which resulted in an improved description of shock structure relative to the Navier-Stokes 
solution - surface boundary conditions were obtained from order-An slip relations. The 
rationale for this approach was that the flow adjacent to the surfaces (within the Prandtl 
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boundary layer) could be well described in the Navier-Stokes level, and Burnett effects 
would only become significant within the shock wave itself. 10 

This is not the case in the problem examined here, as thermal stress will affect the 
pressure field throughout the heat flow domain. However, following the rationale used to 
develop the inner solution, we submit that an order -Kn pressure slip relation can be used 
as a boundary condition, providing that it is corrected for the order-An 2 bulk thermal 
stress effect in the Knudsen layer. Specifically, the quantity / 1C appearing in Eq. (21) 
would be related to the pressure slip at the cold surface via 


fic = 


Kn 


(<pC ~ 4>C, oo) 


Kn—vO 


0C,oo Cp 


(24) 


in which <j > c , oc is the ‘freestream’ pressure adjacent to the cold gas, d' c is the temperature 
gradient at the surface, and c p is a constant. We now assume that <f> c , 0 0 is given by the 
bulk distribution averaged across the Knudsen layer, per the procedure in Eqs. (19) and 
(21). This results in 


^c = r‘(\-Knc p e' c ) + Kn 7 {^KA\ (25) 

The value of c p can be estimated from the results of Loyalka 19 and Sone et al. 12 , who 
calculated the distributions of temperature and density in the Knudsen layer of a heated, 
hard-sphere gas by numerical solution of the linearized Boltzmann equation. The pressure 
at the surface would be obtained from the product of temperature and density at the 
surface, which, to order Kn, yields c p = 0.1913. 


III. NUMERICAL COMPUTATION 
A. Direct Simulation Monte Carlo method 

The 1-D DSMC procedure used to calculate the density, temperature, and pressure 
is based on a code used in an earlier investigation 15 , and follows that developed by Bird. 25 
Because of the highly nonisothermal conditions of the simulations, a nonuniform grid was 
employed in which the length Ax of each cell (or control volume) was set at 0.025-0.1 of 
the local mean free path as estimated from the continuum temperature profile. This is 
equivalent to maintaining a nominally constant number of computational molecules per 
cell. Typically, 20 molecules were assigned per cell, and the time step in the simulations, 
At, was 0.0125-0.05 of the average collision time. Collisions between pairs of molecules in 
a cell were simulated using the no-time-counter approach. 


10 



Calculations proceeded from an initial equilibrium state corresponding to the refer 
ence temperature and pressure of the gas. Sampling of molecular properties was begun 
after a fixed period (typically around 1000 collision times) from the starting point, to elim- 
inate bias in the final results from the initial state. Properties of the bulk gas (density, 
temperature, pressure) were obtained from time averages of molecular properties sampled 
from each cell. The heat flux and normal stress were obtained from recording of the net 
molecular momentum and kinetic energy transport at the boundaries. 

All DSMC calculations were performed on a SUN Ultrasparc I workstation. Presented 
results correspond to averages taken over 2 x 10 7 — 4 x 10 7 times steps which required 
run times of 3-5 days on our machine. The large amount of averaging was required to 
reduce the noise (i.e., natural fluctuations) in the pressure profiles to a point that allowed 
comparison with the theoretical, continuum-based predictions. 

B. Numerical Solution of Continuum Model 

A finite difference method was used to numerically solve Eqs. (13) and (25) (and the 
corresponding boundary condition at the heated surface) for the pressure distribution. The 
domain was discretized onto a uniform mesh, and central differences were used to represent 
the pressure derivatives. The resulting difference equations were solved iteratively for the 
nodal values of 0, using an initial value of r* estimated from Eq. (17). The integral of 
(j> over the domain was then computed using a trapezoid rule, and this value was divided 
into r* to obtain a new estimate of the normal stress. The procedure was then repeated 
until the relative change in t* decreased below 10 Note that the convergence in t to a 
constant implies that the calculated pressure field obeys Eq. (15). The solution procedure 
was stable for all values of Kn and q* presented here, and the calculated values of <f> and 
t* were independent of the initial values used in the iterations. Typically, 200 mesh points 
were used in the calculations. In all cases, presented results are within 0.5% of the values 
obtained when the number of mesh points was doubled. 

IV. RESULTS 

We begin by examining the behavior of the gas that is contained between two solid, 
perfectly-accommodating surfaces. Dimensionless wall temperatures, for the DSMC cal- 
culations, were calculated from Eqs. (8) and (10) for set values of q* and Kn. Three values 
of q* (= 1.0, 1.5, and 2.0) and three of Kn (= 0.05, 0.10, and 0.20) were used in the 
simulations. Listed in Table I are the values of 9cs, @H, an d ^H5i calculated from the 
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slip/continuum relations, corresponding to the q* and Kn combinations. 

The simulations presented here can represent highly nonequilibrium conditions. In the 
case of q* = 1.5 and Kn = 0.2, the ratio 9hs/Qcs ~ which represents the ratio of the hot 
and cold surface temperatures - is around 23. Although such temperature gradients are 
obviously not representative of ‘common’ heat transfer applications, we found it necessary 
to include these cases to create situations in which the effects of thermal stress on pressure 
became clearly apparent. 

Figure 2 presents temperatures, calculated from DSMC simulations, for q* = 1.5 and 
Kn = 0.05, 0.1 and 0.2. Shown also is the continuum prediction from Eq. (7). The 
dimensionless continuum prediction is independent of Kn - since it is based on the heat 
flux in the gas and not on the surface temperatures. The results are entirely consistent 
with previous DSMC investigations of 1-D heat transfer. 17 ’ 18 The DSMC and continuum 
temperatures are in good agreement in the bulk gas - which supports the accuracy of 
the temperature slip conditions. The breakdown of the continuum model is also clearly 
evident within the Knudsen layer adjacent to the hot wall - which, as expected, increases 
approximately linearly in thickness with increasing Kn. 

Comparison of the DSMC and Burnett predictions of pressure distribution appear in 
Figs. 3-5. The quantity <p/r * = P/r, where r is the DSMC-calculated normal stress, is 
plotted vs. £ for q* = 1.0 (Fig. 3). 1.5 (Fig. 4), and 2.0 (Fig. 5). The Knudsen number Kn 
is a parameter in each plot. We present (p/r* , rather than 0, in order to separate the curves 
and to illustrate the disparity between the pressure and the normal stress. Two theoretical 
predictions of P/r are shown for each set of DSMC results. The first corresponds to the 
numerical solution of Eq. (13) with boundary conditions given by Eq. (25) and temper- 
atures predicted from Eq. (7). The second represents the outer approximation given by 
Eq. (14), in which t* is predicted by Eq. (17). This second curve has been adjusted by a 
constant to match the numerical solution results at £ = 0.5. 

As is evident from the DSMC results, the pressure profiles show distinct Knudsen 
layers at both the cooled and heated surfaces, with pressure decreased and increased rela- 
tive to the freestream values at the cold and hot surfaces, respectively. The pressure drop 
at the cold surface can be considerable for the conditions examined here - amounting to 
around an 8% decrease for Kn = 0.2 and q* = 1.5. 

The numerical solution of the Burnett equation, with the corrected slip boundary 
conditions, is seen to capture the essential features of the DSMC pressure distribution. 
In particular, the solution provides a reasonable description of the width and form of the 
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Knudsen layers and the pressure distribution outside the layers. The difference between 
the theoretical and DSMC results is greatest at the edge of the cold-surface Knudsen 
layers, and, as expected, this difference increases with increasing Kn and q* . The same 
trends are seen in the differences between the DSMC pressures at the surface and those 
calculated by Eq. (25), yet this difference is no larger than 10%, on a relative basis, 
for the conditions examined here. We also examined solutions to Eq. (13) in which the 
thermal stress correction was removed from Eq. (25), and these results were in considerably 
less agreement with the DSMC values. Furthermore, solutions to Eq. (13) that used 
boundary conditions derived from the DSMC pressures at the surface (i.e., ‘exact’ boundary 
conditions) did not result in a significantly better fit with the DSMC pressure distributions 
than those generated with Eq. (25). 

A comparison of the DSMC, numerical Burnett solution, and composite inner asymp- 
totic approximation of Eq. (23) is presented in Fig. 6, for the conditions given in Fig. 4 
(i.e., q* = 1.5). We use a value of r* in Eq. (23) corresponding to that predicted from 
Eq. (17) and boundary conditions at £ = 0 are obtained from Eq. (25). Note that Eq. (15) 
would not hold for Eq. (23), since the solution cannot match both boundary conditions. 
Results are shown for £ = 0 to 0.2. The inner and full-Burnett solutions are practically 
indistinguishable for Kn = 0.05 and 0.1, yet for Kn = 0.2 there is a substantial diver- 
gence between the two results in the freestream. This latter difference indicates that ‘bulk’ 
conditions are not attained for the given conditions, i.e., the Knudsen layers affect the 
pressure distribution throughout the heat flow domain. The same conclusion can be in- 
ferred from comparison of the outer approximations and Burnett predictions in Figs. 3 and 
4 for Kn = 0.2. On the other hand, the outer distribution in pressure is evident in the 
results for smaller Kn (i.e., Fig. 5). 

As can be theoretically inferred from Eq. (23), as well as by inspection of Figs. 3 and 5, 
decreasing Kn while maintaining Kn q* constant results in progressively narrower Knudsen 
layers, for which the effects of thermal stress in the bulk gas can be better observed. It 
would be impractical to perform DSMC calculation for values of q* much larger than the 
maximum of 2.0 used here. We can, however, alter the DSMC boundary conditions at 
£ = 1 to model molecular conditions in the bulk gas, as opposed to molecular reflection 
from a solid surface. This approach will, in principle, eliminate the Knudsen layer effects 
at £ = 1, and thus better isolate the effects of thermal stress. 

Following this rationale, we performed additional DSMC calculations on the ‘open’ 
boundary system. To maintain zero net mass flux at the £ = 1 boundary, a new molecule 
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was released into the system every time a molecule left the system by crossing £ = 1. The 
velocity components of the released molecule were randomly sampled from a Chapman- 
Enskog distribution function, with probability given by 


p{u\s •) = 


1 - Du * 


u 2 + s* 2 


-d: 


(26) 


In the above, 


u 


(27) 


is the dimensionless axial thermal velocity of the molecule, s* = y/v * 2 + ie* 2 is the dimen- 
sionless radial molecular velocity, and 


15-v/tt A dT 
16 Th dx 


15\/7r Knq* 

16 n 1 / 2 

°H 


(28) 


Note that the constant D is chosen to simulate a heat flux of q* entering through the hot 
boundary, which is maintained at a temperature of 9h- An acceptance-rejection procedure 
was used to sample velocities from Eq. (26). Aside from the changed distribution function, 
the computational procedure remained the same as before. Calculations were performed 
for Kn = 0.1 (i.e., the system outer boundary is located at 10A re /), and q* was set to 1.0, 
1.5, and 2.0. The value of 9h used in Eq. (26) is the same as that listed in Table I for the 
corresponding value of q* . For q* = 2.0, the value of D is 0.228 - which is near the limit 
of validity of the Chapman-Enskog distribution function. 25 

DSMC results of temperature for the open-boundary system, and the continuum 
predictions given by Eq. (7), appear in Fig. 7. Knudsen layers are no longer observed 
at £ = 1, which is indicated by the close correspondence of the DSMC and continuum 
predictions throughout the system. Pressure distribution results are plotted in Fig. 8. 
Theoretical predictions again correspond to the full solution of Eq. (13) and the outer 
approximation of Eq. (14). The boundary condition at £ = 1 for the numerical solution 
has now been changed to remove the slip effect (i.e., c p — 0 in Eq. (25) for £ = 1). As was 
the case with the temperature profiles, the DSMC pressure distributions do not show the 
presence of Knudsen layer effects at £ = 1. For q* < 1.5, the DSMC pressure distribution 
in the bulk gas is seen to correspond closely to the prediction of Eq. (14) - which in turn 
matches the distribution predicted by the full solution of Eq. (13). These results indicate 
that the bulk-gas pressure distribution is a direct effect of temperature gradients within 
the gas - and is not the result of nonequilibrium at the surfaces. 
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A final comparison of theory and experiment can be obtained from the dimensionless 
normal stress. The simplified model of Eq. (17) indicates that r* should be a function 
primarily of Knq*. Accordingly, we plot in Fig. 9 the DSMC values of r* vs. Knq* for the 
seven different combinations of Kn and q* that were used in the closed system calculations 
(results of Figs. 3-5). Theoretical results correspond to the derived eigenvalues of Eq. (13) 
and to the approximation given by Eq. (17). 

The first point to make is that the predictions of r* from full solution of the Burnett 
differential equation are closely equivalent to those obtained from the outer approximation 
of Eq. (17). Evidently, the decrease in pressure at the cold surface is compensated by 
the increase at the hot, so that the Knudsen layers have a small effect on the averaged 
pressure in the gas. Secondly, the primary dependence of r* on Knq* is supported in 
the DSMC results at Knq* = 0.1 and 0.2 - which each correspond to two combinations 
of Kn and q*. As observed, the results are nearly identical at these points. Finally, the 
theoretical predictions are in excellent agreement with the DSMC results for Knq* < 0.15, 
beyond which the theory overpredicts the decrease in r* . As can be seen from the results, 
the relative decrease in normal stress on the surfaces is quite small, i.e., r* = 0.975 for 
q* = 2.0 and Kn = 0.2, or a 2.5% decrease in ‘measured’ pressure at the surface. We 
should emphasize, however, that this decrease is still significantly larger than the numerical 
precision of the DSMC simulations. 

V. SUMMARY 

This investigation has compared Burnett equation predictions of pressure distribution 
and normal stress with DSMC results for 1-D heat transfer in a stationary, hard-sphere gas. 
For Knq* less than around 0.1 - 0.15, our results show that the theoretical pressure profiles 
can reasonably match those obtained from DSMC. The pressure distribution within the 
Knudsen layer is consistent with an order- Ah asymptotic solution of the Burnett equation. 
Likewise, a simplified outer solution to the Burnett equation, which neglects Knudsen layer 
effects at the surfaces, gives an accurate prediction of the reduction in normal stress due 
to thermal stress, for the same constraint on Knq*. 

We conclude by emphasizing that our findings do not imply that the Burnett equations 
and the derived pressure boundary conditions will remain valid when applied to more 
complicated heat and momentum transfer configurations - such as those involving gas 
flows and velocity gradients. Indeed, the stationary, 1-D system examined here probably 
represents the most simple configuration in which thermal stress effects could be isolated. 
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In this sense, a fatal flaw in the Burnett equations would have been revealed by a lack 
of correspondence between the experimental (i.e., DSMC) observations and theoretical 
predictions of pressure and stress in this system. The results of this investigation, however, 
prove otherwise, and hopefully will encourage others to examine the veracity and usefulness 
of the Burnett equations in practical rarefied gas dynamics applications. 
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q * 

Kn 

®cs 

0 c 

Oh 

&hs 

1.0 

0.10 

0.4451 

0.5469 

1.5364 

1.7791 

1.0 

0.20 

0.3753 

0.5469 

1.5364 

2.1129 

1.5 

0.05 

0.2986 

0.3615 

1.8260 

2.0151 

1.5 

0.10 

0.2543 

0.3615 

1.8260 

2.2479 

1.5 

0.20 

0.1961 

0.3615 

1.8260 

2.9235 

2.0 

0.05 

0.1528 

0.2092 

2.1241 

2.4029 

2.0 

0.10 

0.1203 

0.2092 

2.1241 

2.7659 


TABLE I. Dimensionless surface and gas temperatures. 


Figure Captions 


1. Problem schematic 

2. DSMC and continuum dimensionless temperature distributions vs. dimensionless po- 
sition, q* = 1.0. 

3. DSMC and continuum pressure distributions, vs. dimensionless position, q* = 1.0. 

4. DSMC and continuum pressure distributions, vs. dimensionless position, q* = 1.5. 

5. DSMC and continuum pressure distributions, vs. dimensionless position, q* = 2.0. 

6. DSMC, exact Burnett, and asymptotic approximation pressure distributions, vs. di- 

mensionless position, q* = 1.5. 

7. DSMC and continuum temperature distributions, open system. 

8. DSMC and continuum pressure distributions, open system. 

9. DSMC and continuum predictions of dimensionless normal stress, vs. Knq*. 
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